rm(list=ls())
library(locfit)
setwd("/Users/innoshuadmin/Box Sync/DRD")
source("./functions/DRDfunctions.R")
## load in robust bandwidth selection functions from CCT(2014)## 
source("./functions/rdbwselect.R")
source("./functions/functions.R")
source("./simulations/simuls_brooklyn.R") 

numn=4
nummodel=1
numsm=6
numsim=500
MRDstat=matrix(NA,nrow=numn*nummodel,ncol=numsm)
DRDstat=matrix(NA,nrow=numn*nummodel,ncol=numsm)
DRDstat2=matrix(NA,nrow=numn*nummodel,ncol=numsm)
MRDbw=matrix(NA,nrow=numn*nummodel,ncol=numsm)
DRDbw=matrix(NA,nrow=numn*nummodel,ncol=numsm)
DRDbw2=matrix(NA,nrow=numn*nummodel,ncol=numsm)
size1=matrix(NA,nrow=numn*nummodel,ncol=numsm)
size2=matrix(NA,nrow=numn*nummodel,ncol=numsm)
power1=matrix(NA,nrow=numn*nummodel,ncol=numsm)
power2=matrix(NA,nrow=numn*nummodel,ncol=numsm)
validsim=matrix(NA,nrow=numn*nummodel,ncol=numsm)

for (nind in 1:numn){
	n=1000*2^(nind-1)
	#for (model in 1:nummodel){
		model=1
		for (smind in 1:numsm){
			print(c(nind,model,smind))
			if (smind<=3){
				smooth=3.5+0.5*smind
				bwsel="CCT"
			} else {
				smooth=3+0.5*(smind-3)
				bwsel="IK"				
			}
			set.seed(123)
			MRDbw.sim=rep(NA,numsim)
			MRDstat.sim=rep(NA,numsim)
			DRDbw.sim=rep(NA,numsim)
			DRDstat.sim=rep(NA,numsim)
			DRDbw2.sim=rep(NA,numsim)
			DRDstat2.sim=rep(NA,numsim)
			size1.sim=rep(NA,numsim)
			size2.sim=rep(NA,numsim)
			power1.sim=rep(NA,numsim)
			power2.sim=rep(NA,numsim)
				for (s in 1:numsim){
					simulation_brooklyn
				results=simulation_brooklyn(model,n,smooth,bwsel=bwsel)
				MRDbw.sim[s]=results$MRDbw
				MRDstat.sim[s]=abs(results$MRDstat)
				DRDbw.sim[s]=results$DRDbw
			    DRDstat.sim[s]=max(abs(na.omit(results$DRDstat)))
				DRDbw2.sim[s]=results$DRDbw2
			    DRDstat2.sim[s]=max(abs(na.omit(results$DRDstat2)))
			    size1.sim[s]=sum(results$rightsettrue*0+1)-sum(results$set1 %in% results$rightsettrue)
			    size2.sim[s]=sum(results$set2*0+1)-sum(results$set2 %in% results$rightsettrue)
			    power1.sim[s]=sum(results$set1 %in% results$rightsettrue)/sum(results$set1*0+1)
			    power2.sim[s]=sum(results$set2 %in% results$rightsettrue)/sum(results$rightsettrue*0+1)			    	
				}
		validsim[(nind-1)*nummodel+model,smind]=sum(na.omit(DRDstat.sim*0+1))
		MRDbw[(nind-1)*nummodel+model,smind]=mean(na.omit(MRDbw.sim))
		MRDstat[(nind-1)*nummodel+model,smind]=mean(na.omit(MRDstat.sim)>qt(0.975,n))
		DRDbw[(nind-1)*nummodel+model,smind]=mean(na.omit(DRDbw.sim))
		DRDstat[(nind-1)*nummodel+model,smind]=mean(na.omit(DRDstat.sim)>1.3581)
		DRDbw2[(nind-1)*nummodel+model,smind]=mean(na.omit(DRDbw2.sim))
		DRDstat2[(nind-1)*nummodel+model,smind]=mean(na.omit(DRDstat2.sim)>1.3581)
		size1[(nind-1)*nummodel+model,smind]=mean(na.omit(size1.sim)>0)
		size2[(nind-1)*nummodel+model,smind]=mean(na.omit(size2.sim)>0)
		power1[(nind-1)*nummodel+model,smind]=mean(na.omit(power1.sim))
		power2[(nind-1)*nummodel+model,smind]=mean(na.omit(power2.sim))
		save.image(file=paste("simu_model_",model,".RData",sep=""))
		}
#	}
print(MRDstat)
print(DRDstat)
print(DRDstat2)
print(MRDbw)
print(DRDbw2)
print(size1)
print(size2)
print(power1)
print(power2)
print(validsim)	
}


  